MODULE SRFWDIFSAD_MOD
CONTAINS
SUBROUTINE SRFWDIFSAD(KIDIA,KFDIA,KLON,KLEVS,&
 & PWSAM1M5 ,PCFW5  ,PRHSW5   ,PCDZ5 ,&
 & PWSADIF5 ,LDLAND ,LDALLAYS ,YDSOIL,&
 & PWSAM1M  ,PCFW   ,PRHSW    ,PCDZ  ,&
 & PWSADIF  )

USE PARKIND1  , ONLY : JPIM, JPRB
USE YOMHOOK   , ONLY : LHOOK, DR_HOOK, JPHOOK
USE YOS_SOIL  , ONLY : TSOIL

#ifdef DOC
! (C) Copyright 2012- ECMWF.
!
! This software is licensed under the terms of the Apache Licence Version 2.0
! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
! In applying this licence, ECMWF does not waive the privileges and immunities
! granted to it by virtue of its status as an intergovernmental organisation
! nor does it submit to any jurisdiction.

!**** *SRFWDIFSAD* -  DOES THE IMPLICIT CALCULATION FOR SOIL MOISTURE
!                     (Adjoint)

!     PURPOSE.
!     --------
!          SOLVE TRIDIAGONAL SYSTEM OF EQUATIONS FOR SOIL MOISTURE.

!**   INTERFACE.
!     ----------
!          *SRFWDIFSAD* IS CALLED FROM *SRFTSAD,SRFISAD*

!     PARAMETER   DESCRIPTION                                    UNITS
!     ---------   -----------                                    -----

!     INPUT PARAMETERS (INTEGER):
!     *KIDIA*      START POINT
!     *KFDIA*      END POINT
!     *KLON*       NUMBER OF GRID POINTS PER PACKET
!     *KLEVS*      NUMBER OF SURFACE LAYERS

!     INPUT PARAMETERS (LOGICAL):
!     *LDLAND*     LAND/SEA MASK (TRUE/FALSE)
!     *LDALLAYS*   TRUE - COMPUTATION FOR ALL LAYERS , FALSE - FOR TOP ONLY

!     INPUT PARAMETERS AT T-1 (REAL):
!  Trajectory  Perturbation  Description                                Unit
!  PWSAM1M5    PWSAM1M       PSI AT TIME LEVEL T
!  PCFW5       PCFW          MODIFIED DIFFUSIVITIES (LAMBDA-STAR IN DOCUM.)
!                            (INDEX JK REFERS TO LEVEL JK+1/2)
!  PRHSW5      PRHSW         RIGHT-HAND SIDE OF EQUATIONS
!                            TREATMENT OF EVAPORATION FLUX
!  PCDZ5       PCDZ          C*DZ  (PROFILE OF LAYER DEPTH FOR SOIL WATER)

!     OUTPUT PARAMETERS (REAL):
!  Trajectory  Perturbation  Description                                Unit
!  PWSADIF5    PWSADIF       PSI-STAR  DIVIDED BY ALFA

!     METHOD.
!     -------
!     *LU*-DECOMPOSITION AND BACK SUBSTITUTION IN ONE DOWNWARD SCAN
!     AND ONE UPWARD SCAN.

!     EXTERNALS.
!     ----------
!          NONE.

!     REFERENCE.
!     ----------
!          SEE SOIL PROCESSES' PART OF THE MODEL'S DOCUMENTATION FOR
!     DETAILS ABOUT THE MATHEMATICS OF THIS ROUTINE.

!     Original
!     --------
!     M. Janiskova              E.C.M.W.F.     02-04-2012

!     Modifications
!     -------------

!     ------------------------------------------------------------------
#endif

IMPLICIT NONE

! Declaration of arguments

INTEGER(KIND=JPIM), INTENT(IN)    :: KIDIA
INTEGER(KIND=JPIM), INTENT(IN)    :: KFDIA
INTEGER(KIND=JPIM), INTENT(IN)    :: KLON
INTEGER(KIND=JPIM), INTENT(IN)    :: KLEVS

LOGICAL,            INTENT(IN)    :: LDLAND(:)
LOGICAL,            INTENT(IN)    :: LDALLAYS

TYPE(TSOIL),        INTENT(IN)    :: YDSOIL

REAL(KIND=JPRB),    INTENT(IN)    :: PWSAM1M5(:,:)
REAL(KIND=JPRB),    INTENT(IN)    :: PCFW5(:,:)
REAL(KIND=JPRB),    INTENT(IN)    :: PRHSW5(:,:)
REAL(KIND=JPRB),    INTENT(IN)    :: PCDZ5(:,:)
REAL(KIND=JPRB),    INTENT(INOUT) :: PWSADIF5(:,:)

REAL(KIND=JPRB),    INTENT(INOUT):: PWSAM1M(:,:)
REAL(KIND=JPRB),    INTENT(INOUT) :: PCFW(:,:)
REAL(KIND=JPRB),    INTENT(INOUT) :: PRHSW(:,:)
REAL(KIND=JPRB),    INTENT(INOUT) :: PCDZ(:,:)
REAL(KIND=JPRB),    INTENT(INOUT) :: PWSADIF(:,:)

!*         0.2    DECLARATION OF LOCAL VARIABLES.
!                 ----------- -- ----- ----------

REAL(KIND=JPRB) :: ZTCOE5(KLON,KLEVS), ZEBSW5(KLON,KLEVS)
REAL(KIND=JPRB) :: ZQDP5(KLON,KLEVS), ZDISC5(KLON,KLEVS)
REAL(KIND=JPRB) :: ZWSADIFA5(KLON,KLEVS),ZWSADIFB5(KLON,KLEVS)
REAL(KIND=JPRB) :: ZFAC5(KLON,KLEVS)

REAL(KIND=JPRB) :: ZTCOE(KLON), ZEBSW(KLON,KLEVS)
REAL(KIND=JPRB) :: ZDISC, ZFAC, ZQDP, ZTPFAC2

INTEGER(KIND=JPIM) :: ILEVM1, JK, JL
REAL(KIND=JPHOOK) :: ZHOOK_HANDLE

!     ------------------------------------------------------------------

!*         1.    SET UP SOME CONSTANTS.
!                --- -- ---- ----------
IF (LHOOK) CALL DR_HOOK('SRFWDIFSAD_MOD:SRFWDIFSAD',0,ZHOOK_HANDLE)
ASSOCIATE(RSIMP=>YDSOIL%RSIMP)

ZTPFAC2=1.0_JPRB/RSIMP
ILEVM1=KLEVS-1

!*         1.1     SETTING OF RIGHT HAND SIDES.

DO JK=1,KLEVS
  DO JL=KIDIA,KFDIA
    IF (LDLAND(JL)) THEN
      PWSADIF5(JL,JK) = ZTPFAC2*PWSAM1M5(JL,JK)+PRHSW5(JL,JK)
      ZWSADIFA5(JL,JK) = PWSADIF5(JL,JK)
    ENDIF
  ENDDO
ENDDO

!*         1.2     TOP LAYER ELIMINATION.

DO JL=KIDIA,KFDIA
  IF (LDLAND(JL)) THEN
    ZTCOE5(JL,1) = PCFW5(JL,1)
    ZQDP5(JL,1) = 1.0_JPRB/PCDZ5(JL,1)
    ZDISC5(JL,1) = 1.0_JPRB/(1.0_JPRB+PCFW5(JL,1)*ZQDP5(JL,1))
    ZEBSW5(JL,1) = ZDISC5(JL,1)*(PCFW5(JL,1)*ZQDP5(JL,1))
    PWSADIF5(JL,1) = ZDISC5(JL,1)*ZWSADIFA5(JL,1)
    ZWSADIFB5(JL,1) = PWSADIF5(JL,1)
  ENDIF
ENDDO


IF (LDALLAYS) THEN
!*         1.3     ELIMINATION FOR MIDDLE LAYERS.

  DO JK=2,ILEVM1
    DO JL=KIDIA,KFDIA
      IF (LDLAND(JL)) THEN
        ZQDP5(JL,JK) = 1.0_JPRB/PCDZ5(JL,JK)
        ZFAC5(JL,JK) = ZTCOE5(JL,JK-1)*ZQDP5(JL,JK)
        ZTCOE5(JL,JK) = PCFW5(JL,JK)
        ZDISC5(JL,JK) = &
         &  1.0_JPRB/(1.0_JPRB+ZFAC5(JL,JK)*(1.0_JPRB-ZEBSW5(JL,JK-1)) &
         & +PCFW5(JL,JK)*ZQDP5(JL,JK))
        ZEBSW5(JL,JK) = ZDISC5(JL,JK)*(PCFW5(JL,JK)*ZQDP5(JL,JK))
        PWSADIF5(JL,JK) = ZDISC5(JL,JK) &
         & *(ZWSADIFA5(JL,JK)+ZFAC5(JL,JK)*ZWSADIFB5(JL,JK-1))
        ZWSADIFB5(JL,JK) = PWSADIF5(JL,JK)
      ENDIF
    ENDDO
  ENDDO

!*         1.4     BOTTOM LAYER ELIMINATION.

  DO JL=KIDIA,KFDIA
    IF (LDLAND(JL)) THEN
      ZQDP5(JL,KLEVS) = 1.0_JPRB/PCDZ5(JL,KLEVS)
      ZFAC5(JL,KLEVS) = ZTCOE5(JL,KLEVS-1)*ZQDP5(JL,KLEVS)
      ZDISC5(JL,KLEVS) = &
       &  1.0_JPRB/(1.0_JPRB+ZFAC5(JL,KLEVS)*(1.0_JPRB-ZEBSW5(JL,KLEVS-1)) &
       & +PCFW5(JL,KLEVS)*ZQDP5(JL,KLEVS))
      PWSADIF5(JL,KLEVS) = ZDISC5(JL,KLEVS) &
       & *(ZWSADIFA5(JL,KLEVS)+ZFAC5(JL,KLEVS)*ZWSADIFB5(JL,KLEVS-1))
      ZWSADIFB5(JL,KLEVS) = PWSADIF5(JL,KLEVS)
    ENDIF
  ENDDO

!*         1.5     BACK-SUBSTITUTION.

  DO JK=ILEVM1,1,-1
    DO JL=KIDIA,KFDIA
      IF (LDLAND(JL)) THEN
        PWSADIF5(JL,JK) = PWSADIF5(JL,JK)+ZEBSW5(JL,JK)*PWSADIF5(JL,JK+1)
      ENDIF
    ENDDO
  ENDDO
ENDIF

!          0.  ADJOINT CALCULATIONS
!              --------------------

!* Set local variables to zero

ZTCOE(:) = 0.0_JPRB
ZEBSW(:,:) = 0.0_JPRB

IF (LDALLAYS) THEN

!*         0.1.5     BACK-SUBSTITUTION.

  DO JK=1,ILEVM1
    DO JL=KIDIA,KFDIA
      IF (LDLAND(JL)) THEN
        PWSADIF(JL,JK+1) = PWSADIF(JL,JK+1)+ZEBSW5(JL,JK)*PWSADIF(JL,JK)
        ZEBSW(JL,JK) = ZEBSW(JL,JK)+PWSADIF5(JL,JK+1)*PWSADIF(JL,JK)
      ENDIF
    ENDDO
  ENDDO

!*         0.1.4     BOTTOM LAYER ELIMINATION.

  DO JL=KIDIA,KFDIA
    IF (LDLAND(JL)) THEN
      ZQDP = 0.0_JPRB
      ZFAC = 0.0_JPRB
      ZDISC = 0.0_JPRB

      ZDISC = ZDISC+(ZWSADIFA5(JL,KLEVS)+ZFAC5(JL,KLEVS)*ZWSADIFB5(JL,KLEVS-1)) &
       & *PWSADIF(JL,KLEVS)
      PWSADIF(JL,KLEVS-1) = PWSADIF(JL,KLEVS-1) &
       & +ZDISC5(JL,KLEVS)*ZFAC5(JL,KLEVS)*PWSADIF(JL,KLEVS)
      ZFAC = ZFAC+ZDISC5(JL,KLEVS)*ZWSADIFB5(JL,KLEVS-1)*PWSADIF(JL,KLEVS)
      PWSADIF(JL,KLEVS) = ZDISC5(JL,KLEVS)*PWSADIF(JL,KLEVS)

      ZFAC = ZFAC-(1.0_JPRB-ZEBSW5(JL,KLEVS-1)) &
       & *ZDISC5(JL,KLEVS)*ZDISC5(JL,KLEVS)*ZDISC
      ZEBSW(JL,KLEVS-1) = ZEBSW(JL,KLEVS-1) &
       & +ZFAC5(JL,KLEVS)*ZDISC5(JL,KLEVS)*ZDISC5(JL,KLEVS)*ZDISC
      ZQDP = ZQDP-PCFW5(JL,KLEVS)*ZDISC5(JL,KLEVS)*ZDISC5(JL,KLEVS)*ZDISC
      PCFW(JL,KLEVS) = PCFW(JL,KLEVS) &
       & -ZQDP5(JL,KLEVS)*ZDISC5(JL,KLEVS)*ZDISC5(JL,KLEVS)*ZDISC

      ZQDP = ZQDP+ZTCOE5(JL,KLEVS-1)*ZFAC
      ZTCOE(JL) = ZTCOE(JL)+ZQDP5(JL,KLEVS)*ZFAC
      PCDZ(JL,KLEVS) = PCDZ(JL,KLEVS)-ZQDP5(JL,KLEVS)*ZQDP5(JL,KLEVS)*ZQDP
    ENDIF
  ENDDO
      
!*         0.1.3     ELIMINATION FOR MIDDLE LAYERS.

  DO JK=ILEVM1,2,-1
    DO JL=KIDIA,KFDIA
      IF (LDLAND(JL)) THEN
        ZQDP = 0.0_JPRB
        ZFAC = 0.0_JPRB
        ZDISC = 0.0_JPRB

        ZDISC = ZDISC+(ZWSADIFA5(JL,JK)+ZFAC5(JL,JK)*ZWSADIFB5(JL,JK-1)) &
         & *PWSADIF(JL,JK)
        PWSADIF(JL,JK-1) = PWSADIF(JL,JK-1) &
         & +ZDISC5(JL,JK)*ZFAC5(JL,JK)*PWSADIF(JL,JK)
        ZFAC = ZFAC+ZDISC5(JL,JK)*ZWSADIFB5(JL,JK-1)*PWSADIF(JL,JK)
        PWSADIF(JL,JK) = ZDISC5(JL,JK)*PWSADIF(JL,JK)

        ZDISC = ZDISC+(PCFW5(JL,JK)*ZQDP5(JL,JK))*ZEBSW(JL,JK)
        ZQDP = ZQDP+ZDISC5(JL,JK)*PCFW5(JL,JK)*ZEBSW(JL,JK)
        PCFW(JL,JK) = PCFW(JL,JK)+ZDISC5(JL,JK)*ZQDP5(JL,JK)*ZEBSW(JL,JK)
        ZEBSW(JL,JK) = 0.0_JPRB

        ZFAC = ZFAC-(1.0_JPRB-ZEBSW5(JL,JK-1)) &
         & *ZDISC5(JL,JK)*ZDISC5(JL,JK)*ZDISC
        ZEBSW(JL,JK-1) = ZEBSW(JL,JK-1) &
         & +ZFAC5(JL,JK)*ZDISC5(JL,JK)*ZDISC5(JL,JK)*ZDISC
        ZQDP = ZQDP-PCFW5(JL,JK)*ZDISC5(JL,JK)*ZDISC5(JL,JK)*ZDISC
        PCFW(JL,JK) = PCFW(JL,JK)-ZQDP5(JL,JK)*ZDISC5(JL,JK)*ZDISC5(JL,JK)*ZDISC

        PCFW(JL,JK) = PCFW(JL,JK)+ZTCOE(JL)
        ZTCOE(JL) = 0.0_JPRB
        ZQDP = ZQDP+ZTCOE5(JL,JK-1)*ZFAC
        ZTCOE(JL) = ZTCOE(JL)+ZQDP5(JL,JK)*ZFAC
        PCDZ(JL,JK) = PCDZ(JL,JK)-ZQDP5(JL,JK)*ZQDP5(JL,JK)*ZQDP
      ENDIF
    ENDDO
  ENDDO

ENDIF

!*         0.1.2     TOP LAYER ELIMINATION.

DO JL=KIDIA,KFDIA
  IF (LDLAND(JL)) THEN
    ZQDP = 0.0_JPRB
    ZDISC = 0.0_JPRB

    ZDISC = ZDISC+ZWSADIFA5(JL,1)*PWSADIF(JL,1)
    PWSADIF(JL,1) = ZDISC5(JL,1)*PWSADIF(JL,1)

    ZDISC = ZDISC+PCFW5(JL,1)*ZQDP5(JL,1)*ZEBSW(JL,1)
    ZQDP = ZQDP+ZDISC5(JL,1)*PCFW5(JL,1)*ZEBSW(JL,1)
    PCFW(JL,1) = PCFW(JL,1)+ZDISC5(JL,1)*ZQDP5(JL,1)*ZEBSW(JL,1)
    ZEBSW(JL,1) = 0.0_JPRB

    ZQDP = ZQDP-PCFW5(JL,1)*ZDISC5(JL,1)*ZDISC5(JL,1)*ZDISC
    PCFW(JL,1) = PCFW(JL,1)-ZQDP5(JL,1)*ZDISC5(JL,1)*ZDISC5(JL,1)*ZDISC
    PCDZ(JL,1) = PCDZ(JL,1)-ZQDP5(JL,1)*ZQDP5(JL,1)*ZQDP
    PCFW(JL,1) = PCFW(JL,1)+ZTCOE(JL)
    ZTCOE(JL) = 0.0_JPRB
  ENDIF
ENDDO

!*         0.1.1     SETTING OF RIGHT HAND SIDES.

DO JK=KLEVS,1,-1
  DO JL=KIDIA,KFDIA
    IF (LDLAND(JL)) THEN
      PWSAM1M(JL,JK) = PWSAM1M(JL,JK)+ZTPFAC2*PWSADIF(JL,JK)
      PRHSW(JL,JK) = PRHSW(JL,JK)+PWSADIF(JL,JK)
      PWSADIF(JL,JK) = 0.0_JPRB
    ENDIF
  ENDDO
ENDDO

END ASSOCIATE
IF (LHOOK) CALL DR_HOOK('SRFWDIFSAD_MOD:SRFWDIFSAD',1,ZHOOK_HANDLE)

END SUBROUTINE SRFWDIFSAD
END MODULE SRFWDIFSAD_MOD
